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Abstract 

We discuss two important techniques, series expansion and Monte Carlo 
simulation, for random sequential adsorption study. Random sequential ad- 
sorption is an idealization for surface deposition where the time scale of par- 
ticle relaxation is much longer than the time scale of deposition. Particles are 
represented as extended objects which are adsorbed to a continuum surface or 
lattice sites. Once landed on the surface, the particles stick to the surface. We 
review in some details various methods of computing the coverage 9{t) and 
present some of the recent and new results in random sequential adsorption. 

1 Random sequential adsorption problem 

The random sequential adsorption (RSA) problem arises from a number of situations 
in physics, chemistry, and biology. One of the early RSA problem comes from 
polymer chemistry. Flory fij was interested in the dimerization of polymer chain. 
The adjacent monomers on the side of the chain form dimers one at a time until no 
pair of monomers is available. He found the final reacted monomers based purely 
on geometric argument. It turns out that this problem is equivalent to a random 
sequential adsorption of dimers on a one-dimensional lattice. Consider the one- 
dimensional lattice, with all sites empty initially. At time t > 0, dimers, each of 
which occupies two consecutive sites, are dropped on the lattice at a rate k per 
second equally likely at any locations (see Fig. |T](a)). As long as the two sites are 
empty, a dimer lands on the lattice. If any one of the two sites is already occupied, 
the deposition is rejected. This is the basic RSA model. 
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Figure 1: Some typical random sequential adsorption configurations, (a) Dimer de- 
position on a one- dimensional lattice; (b) monomer with nearest neighbor exclusion; 
(c) the jamming state of hard discs on continuum; (d) self-avoiding random walks 
of length N = 5 on square lattice. 



RSA of many different geometric objects has been studied. The deposition of 
unit line segments on a one-dimensional continuum is called car-parking problem 
0. The deposition of discs on two-dimensional continuum || is relevant to the 
adsorption of colloids or protein molecules on glass surface j|, Lattice models 
of various geometries are studied. In Fig. [I] we show some typical configurations of 
RSA models. 

The very basic question of RSA is the time dependence of the coverage, 9(t), 
and the jamming coverage, 6(oo). The approach to the jamming coverage is also 
of considerable interest. It turns out that the approach to the jamming coverage 
is to some extent universal. For lattice models, the approach to jamming is always 
exponential due to the discreteness of the problems || , while for continuum problems 
the asymptotic dependences exhibit power laws which depend on the symmetry of 
the deposited objects [f7j |9j |10 . 
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A comprehensive review on RSA and cooperative sequential adsorptions is given 
by Evans [11]]. Bartelt and Privman jl2| reviewed, among other things, mean-field 
approximations. Progress on experiment work is given by Ramsden Random 
sequential adsorption has many aspects and is rich with problems. In this article, 
we limit our scope to consider two important approximate techniques for computing 
the coverage of RSA. They are the series expansion and Monte Carlo simulation. 
The emphasis is on methods rather than specific results. Only one-dimensional 
problems Jl4|, [15], [H| [17] [18| and quasi-one-dimensional problems ]19|, ^(], ^TJ have 



exact solutions. Two-dimensional problems are most likely intractable. Because of 
this, the approximate methods like series expansion and Monte Carlo are very useful 
for the study of RSA. We present some of the technical details which are usually 
given only briefly in original research papers. With examples, we show how the 
series can be obtained and analyzed with Pade approximation. For Monte Carlo 
simulation, the emphasis is on efficient algorithms for RSA simulation. 



2 Series expansions 



One of the well-developed methods in studying random sequential adsorption is 
series expansion around time t = 0. The method gains experience from series ex- 
pansion studies of other problems in condensed matter physics Wl ] . Series expansion 
approach was the only successful method to obtain nonclassical results in the study 
of critical phenomena for a while before the advent of renormalization group method 
23| . The methods and technique developed, especially the Pade approximation tech- 
nique p4|, have been also applied successfully to random sequential adsorption. The 



technique offers a systematic approach with controlled approximation. Unlike the 
approach with truncated rate equations |25|, |27j, series expansion can be auto- 
mated on computer relatively easily, and the number of terms obtained can be rather 
high. In many cases, the accuracy is higher than or comparable to most accurate 
Monte Carlo simulation results. 

The series expansion for the rate of adsorption as a function of particle density 
has been introduced by Widom [28|. It was shown that this type of expansion is 



similar to virial expansion for the equation of state. In fact, the first two coefficients 
involving one and two-particle distributions are the same. The difference between 
equilibrium system and RSA starts at third order. Series for hard discs on two- 
dimensional continuum up to third order were first obtained by Schaaf and Talbot 
29| . Series expansions for lattice models and continuum systems have been consid- 
ered by Hoffman |30 and Evans [31]]. Evans presented rate equation approach for 
deriving series. The series in time and in density are related by a transformation, 
so they are equivalent. High-order series with the help of computer were derived by 
Baram and Kutasov ||32|| , Dickman et al [33], Bonnier et al [34|, Baram and Fixman 
|35|, and Gan and Wang p6[. 
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2.1 Rate equations 



The RSA process can be described fully with a master equation of the form 

9P({ ' },f) =£r(M,W({4t), (i) 

01 W'} 

where {a} is a set of state variables which completely specify the state of the system; 
P(- • •) is the probability of such a state; and T is a transition matrix. For notational 
convenience, we shall consider discrete lattice models. Continuous space problem 
can be generalized easily. Then {a} denotes a set of occupation variables cr, located 
at site i. 

Although Eq. ([I]) contains the complete information of the process, it is rather 
difficult to deal with efficiently due to the high dimensionality of the problem. Thus, 
in most of the analytic treatments, reduced quantities are worked with. One of the 
most important such quantities is the marginal probability distribution: 

P(G)= £ P(W), (2) 

{cr} , (Tj=(T? for ieG 

where the summation is over all the possible states such that a given set G of sites 
takes a known set of values. That is, P(G) is the probability that a given set of sites 
having specified values. For the standard RSA problem, it is sufficient to consider 
only a set of connected sites which are unoccupied. 

Our task is to write down a set of differential equations for P(G). In principle, 
we can derive them from the master equation (0). This is not really necessary, and 
we can give the equations based on physical meaning of P(G). Let us take a look 
of the equations associated with a random sequential deposition of dimers on a one- 
dimensional lattice. In this case, we consider the probability P n that a consecutive 
n sites are empty. Assuming translational invariance of the initial conditions, this 
probability does not depend on the specific locations. Due to the simple geometry in 
one dimension, the set G of n connected sites being empty can be characterized by 
a single integer n. Let us consider the change of P n by the deposition process. We 
assume that the n empty sites are at i — 1, 2, . . . , n. Clearly, depositions of dimers 
involving the sites i < or i > n do not change P n , while depositions inside empty 
sites on 1 to n destroy the consecutive empty sequence, the probability of which is 
proportional to P n for each of the n — 1 possible ways of depositions. At the two 
ends, the empty sequence can be destroyed with one site at (or n + 1) outside 
the sequence and the other site at 1 (or n) in the sequence, see Fig. 0. In order 
for this to happen, the site just outside the considered sequence should be empty, 
the probability of which is P n +x- There are two ways to do this. Putting all these 
together, we have, 

-k^^j^ = {n-l)P n (t) + 2P n+1 {t), n = 1,2,3, (3) 
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Figure 2: The one-dimensional dimer deposition. The sites 1 to n are known to be 
empty, but the status of other sites are unknown. 

The proportionality constant k sets the time scale. Without loss of generality, we 
can take it to be unity. This equation can be solved |TB| to give 



P n {t) = exp[-(n - l)t - 2 + 2e 
Note that the coverage is just 9(t) — 1 — Pi{t). 



(4) 



Similar consideration gives so-called rate equations [[37]] for other deposition pro- 
cesses. For example, in deposition of dimers on square lattice, the basic quantity of 
interest is the probability P{G) that the given set G of sites are empty; we do not 
care other sites being occupied or empty. The rate of decrease is proportional to 
the probability that the current configuration G can be destroyed by depositing a 
dimer with at least one site in G. Thus we have 



dP(G) 
dt 



E P (G') 



(5) 



ways of destroying G 



where the summation runs over all possible ways of depositing a dimer at a pair of 
empty sites such that at least one site is in G; G' = G if two of the dimer sites of a 
deposition attempt are within G, or G' is one site more than G so that a deposition 
with one site in G and one site outside G can be carried out. The first two equations 
look like these: 

dP(o) 



dt 



AP(oo) 



rfP(oo) 
df~ 



P(oo) + 2P(ooo) + 4P(oo). 



(6) 
(7) 



There are four ways to destroy a single empty site, provided that the nearest neighbor 
site is also empty. Using the assumption that initial conditions are lattice symmetry 
invariant, we can write them simply as 4P(oo). The second equation is derived 
similarly. 

For a general discussion, we write the rate equations symbolically as 

dP{G) 



dt 



CP{G), 
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where £ is a linear operator defined by 



CP(G)=J2c G ,P(G'). (9) 

G 1 



The ra-th derivative is then 
d n P(G) 



C n P(G), with P(G)\ t=0 = l for all G. (10) 
We assume that all sites are empty at t — 0. 

2.2 Simple counting methods for lattice models 

We can obtain the rules for series expansion starting from the rate equation, Eq. (§). 
We begin with a concrete example of the random sequential adsorption of single 
site occupation with nearest neighbor exclusion on a square lattice. A site can be 
occupied by a particle, as long as the four nearest neighbors are empty. Figure |l](b) 
shows a sample configuration of the nearest neighbor exclusion RSA model. 

The first few rate equations read: 

j o o oo 

U „,OOO x a r->/ 0000 \ 

o ) = P( o )+4P( oo ), (12) 



o 

oo oo ooo ooo 

,oooo N „„,oooo N ^^.ooooo, , ^,oooo N 



- — Pf oo") = 2Pf oo") + 2P("ooo") + 4P("oo"), (13) 

where a pattern with o's denotes a set of sites of that geometric arrangement which 
are empty (while other sites are unspecificed) . We have taken into account the trans- 
lational and rotational symmetry of the problem. The equations are genenerated as 
follows. For each given configuration G, we consider all possible ways of destroying 
the configuration by a deposition of a particle at each one of the empty sites. If 
the four neighbors are already known to be empty, the rate is simply proportional 
to the original probability P{G). If the site is adjacent to sites of unknown status, 
we require those sites to be empty. Thus this second type of process creates new 
configurations G' which have more empty sites than G. 

From this set of rate equations we can get the n-th derivative of the probability 
that the sites in the set G are empty, P^ n \G) = d n P(G)/dt n \ t= o. With an initial 
condition of empty lattice at t — 0, the zeroth derivative P(G, t = 0) equals 1 for 
all configurations G. The most relevant one is the probability that a single site is 
empty, P(o, t), since 9(t) = 1 — P(o, t) is the coverage. Power series expansion is 
obtained from the derivatives evaluated at t — 0: 

oo ±n 

P(o,t) = £pM(o)-. (14) 

n=0 n - 
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For the nearest neighbor exclusion model, we have P(o) = 1. P'{o) = — 1 from 

o o oo 

Eq. (O). From Eq. (O) and (O), we get P"(o) = -P'{°o°) = P{°o°) + AP{ oo ) = 

o o oo 

5. Similarly, we find P"'{o) = -P"{°o°) = P'{°o°) + AP'{ oo°) = -37. Thus, to 

third order, the expansion is 

P^ t ) = l-t + -t 2 --t 3 + 0(t 4 ). (15) 

The process of obtaining the answer can be thought as counting the number 
of patterns in n generations. We note that for each current configuration G, we 
go over the sites of G once, each generated a next generation pattern, which may 
or may not be the same as the parent pattern. The expansion coefficient S(n) = 
(— l)( n+1 )cp +1 P(o, t)/dt n+1 \ t=0 is simply the total number of possible patterns (or 
sequences), where two-dimensional coordinate xq = (0,0) is fixed at the origin, 
while x\ varies in a domain D(xq), and X2 in a union domain D(xq) U D(x\), . . . , 
and x n varies in a domain D(xq) U D(x\) ■ ■ ■ U D(x n -i). The domain D(x) is the 
center site x plus four nearest neighbor sites of x. We write 



S(n)= £ E E 1- (16) 

zigD(ico) x 2 £D(x )UD(x 1 ) x n £D(x )\jD(x 1 )---\jD(x n - 1 ) 

5(0) is defined by —dP(o,t)/dt\ t=0 . By definition, S(n) is always positive. In terms 
of S(n), the expansion for the rate of adsorption is 

, = _f^ = £; s(n) (=£)!. (17) 

This notation is convenient, for Eq. (^) generalizes to other type of lattice models by 
choosing a different domain D(x) and generalizes to random sequential adsorption 
on continuum. For anisotropic objects, x should be understood as position and 
orientation as well. 

The easiest way of implementing this counting method is to use recursive function 
calls. Consider the following C-like pseudo-code: 

RSA(G, n) 
{ 

S(n) += \G\; 

if (n > N max ) return; 

for each (x E G) { 

RSA(GUD(i),n + l); 

} 

} 

In this program, G is a set of empty sites, which can be represented by a list of 
coordinates of the sites; G U D(x) is the union of the set G and the set consisting of 
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x and its four nearest neighbors; \G\ is the cardinality of the set. The variable G is 
local to each function call, while S(n) is global. We assume that S(n) is initialized 
to zero. The recursion stops after N max in depth, which gives us results for S(0) to 
S(N max ). If we evoke the program by RSA(G, 0), the derivatives d n P(G)/dt n \ t=0 = 
(—l) n S(n — 1), 1 < n < N max + 1, are computed. 

Clearly, the above algorithm can be improved. First of all, we need not count 
a configuration if it is the same as the parent configuration. We shall discuss this 
in detail later in Section 2.5. The second improvement is to rewrite the recursive 
calls by nonrecursive procedure. Since the set G always contains the parent set 
where G is derived, it is possible to use only one list representing all the G's (one 
for each level in the recursion), each ends at some point in the list. Of course, the 
bookkeeping will be more complicated. None of the above suggestions will improve 
the speed of calculation substantially. In next section, we discuss perhaps the most 
efficient algorithm [^8j for the series expansions on lattices. 

2.3 More advanced methods 

We note that in the simple counting algorithm, we did not use the symmetry of the 
problem. Two patterns related by translation or rotation symmetry should have the 
same derivatives, but this is not recognized. In order to exploit fully the symmetry of 
the problem, we use the technique of dynamic programming, and use sophisticated 
data structures. The starting point is again the rate equation, Eq. (|J). It turns out 
that the best algorithm is to combine this method with a simple counting algorithm. 

A general function is coded which returns the right-hand side of Eq. (^) when 
the configuration G, or the set G, is given. The set G is represented as a list of 
coordinates constructed in an ordered manner. Note that each term on the right- 
hand side of the equation is just the probability of the set GUD(i) being empty, 
where x runs over the sites of G. Terms which are identical after translational 
and rotational symmetry consideration are collected as one term with associated 
coefficients. Each rate equation is represented by a node together with a list of 
pointers to other nodes. Each node represents a function characterized by the set 
G. The node contains pointers to the derivatives of this node obtained so far, and 
pointers to the "children" of this node and their associated coefficients, which form 
a symbolic representation of the rate equations. The computation of n-th derivative 
uses the rate equations recursively. Since each node is linked to other nodes, the 
computation of the n-th derivative can be considered as expanding a "tree" (with 
arbitrary number of branches) of depth n. The tree structure for the nearest neighbor 
exclusion model is shown in Fig. |3|. Unlike the standard tree, we allow loops linking 
back to earlier generations. 

The traversal or expansion of the tree |39| can be done in a depth-first fashion 
or a breadth-first fashion. Each has a different computational complexity. A simple 
depth-first traversal requires only a small amount of memory of order n. However, 
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Figure 3: The "tree" structure representing Eq. (|TT|) to (|13"D in the RSA series 
expansion of the nomomer with nearest neighbor exclusion. The arrows point from 
parent to child generations; the numbers indicate the expansion coefficients (number 
of equivalent patterns). 

the time complexity is at least exponential, b n , with a large base b. A breadth- 
first algorithm consumes memory exponentially, even after the number of the rate 
equations has been reduced by taking the symmetry of the problem into account. 
The idea of dynamic programming can be incorporated in the breadth-first expan- 
sion where the intermediate results are stored and referred. To achieve the best 
performance, a hybrid of strategies is used to reduce the computational complexity: 

• Each configuration (pattern) is transformed into its canonical representation; 
all configurations related by lattice symmetry are considered as the same con- 
figuration. 

• We use breadth-first expansion to avoid repeated computations involving the 
same configuration. If a configuration has already appeared in earlier expan- 
sion, a pointer reference is made to the old configuration. Each configuration 
is stored in memory only once. However, storing of all the distinct configura- 
tions leads to a very fast growth in memory consumption. For a quick check 
if a configuration is already stored, we use the standard technique of hashing 

• The last few generations in the tree expansion use a simple depth-first traversal 
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to curb the problem of memory explosion. 
• Parallel computation proves to be useful. 

The program is controlled by two parameters D and C . D is the depth of 
breadth-first expansion of the tree. When depth D is reached, we no longer want to 
continue the normal expansion in order to conserve memory. Instead, we consider 
each leaf node afresh as the root of a new tree. The derivatives up to (n — D)th 
order are computed for this leaf node. The expansions of the leaf nodes are done in 
serial, so that the memory resource can be reused. The parameter C controls the 
number of last C generations which should be computed with a simple depth-first 
expansion algorithm. It is a simple recursive counting algorithm, which uses very 
little memory, and can run fast if the depth C is not very large. In this algorithm 
the lattice symmetry is not treated. 

This technique is used to obtain the RSA series for dimer and nearest neigh- 



bor exclusion models on square and honeycomb lattices ||36|| , RSA with diffusional 
relaxation KJ, and Ising relaxation dynamics ||41||. 



2.4 Series for continuum systems 

The results derived for lattice model can be generalized for continuum system. Con- 
sider the deposition of disc of diameter a = 1. One way to obtain the rules for 
series expansion in continuum is to consider the limit of discretized lattice model. 
The RSA of discs can be thought as approximate lattice problem with a single-site 
deposition and exclusion of sites within some distance of the occupied site. The 
analog of Eq. (|T6| ) is |33| , 



S(ri) = J dx\ J dx 2 ■ ■ ■ J dx n , (18) 

xiED(x ) x 2 eD(x )UD(x 1 ) x n eD(x )UV(xi)-UD(x„-i) 

where x% is interpreted as a d- dimensional vector, dxi is a d- dimensional volume 
element, and x £ D(xi) is the set such that \x — Xi\ < 1. 

The integral in Eq. ( Jl8| ) has a rather complicated integration domain. To simplify 
the integral, we introduce the "Mayer function" of a hard disc system, 

fa = ffa ~ Xj) = { o, ' otherwise'.' ~ ' (19) 

Then we can rewrite the integral with integration domain in the whole space and 
with the integrand consisting of terms of products of fif. 

-S(l)= [ f 01 dx 1} (20) 
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5(2) = J f Q1 dx, |[/ 2 + /i 2 (l + K 



dx? 



and in general 



-l) n S(n) = J dxxj dx 2 ---j dz„n E IK 1 + /<*) 



(21) 



(22) 



fc=l \i=o 



The result of expanding the terms leads to the following graph-theoretic descrip- 
tion: (— l) n S(n) is the sum of the contributions from all connected (n + l)-point 
unlabeled graphs; each term is the value of an integral associated with the graph 
times an integer multiplicity factor. The value of a (labeled) graph is the integral 
of a form //•••/ dx\dx 2 • • • dx n T7 fij, called cluster integral. The numerical factor 
is the number of topologically distinct ways of labeling the graph, subject to the 
following constraint: the labels are put down in sequential order from to n, such 
that the current vertex k is connected to a vertex j < k. There is a one-to-one 
correspondence between a graph and an integral. Each vertex i of the graph cor- 
responds to an integration variable Xi (except vertex 0), and each edge of the 
graph corresponds to a function fy. The definitions differ slightly from the normal 
convention [[EJ. The diagrammatic rules give, for the particle density, the following 
expression for the first few terms, 



P(t) 



t + (o-o)- 



o 0-0/ 3! 

-o 
o-b 



(23) 



This diagrammatic rule is similar to the Mayer theory f|!| for equilibrium fluid 
system. The major difference is the prefactor. In Mayer theory, the integer factor 
is simply the number of different ways of labeling a graph, without imposing an 
ordering. It turns out that the parallel is deeper than this, see Table [l]. Just like the 
topological reduction in Mayer theory, by some transformation, Given [33| arrived 
at an expansion involving only the star graphs (irreducible graphs). Let <fi = dp/dt 
be the rate of adsorption, which is proportional to the area available for deposition. 
We have the following virial expansion 



In. 



n=l 



(24) 



where b n is the sum of contributions from all (n+ l)-point star graphs with the RSA 
multiplicity factor. The first few terms are 



ln0=(o-o)p+( o X o )|+(2 



-i:/:mx()!+o(p 4 ) 



(25) 



The diagrammatic rules can also be applied to lattice models when the Mayer func- 
tion and integrals are properly interpreted. In fact, Hoffman [|30| has derived a 
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expansion variable 


fugacity z 


time t 


quantities 


InH 

pressure pP = — — 
p 1 din- 

9 z V dz 


density p 

, , dp 
adsorption rate <p = — 
at 


basic expansion 
graphs 


In - 00 z n 

71 = 1 

a n = sum of all topologically 
distinct, labeled, connected 
n-point graphs. 


00 ±n 
P = £ «n-T 

n ' 

a n = sum of all topologically 
distinct, labeled, connected 
n-point graphs, subject to 
the labeling constraint that 
a vertex i must be connected 
to a vertex with label less 
than i. 


virial expansion 

topological 
reduction 


ln^ = ln^ = ^6 n ^- 

b n = sum of all topologically 
distinct, labeled, doubly 
connected (irreducible) 
(n + l)-point graphs. 

00 K p n+1 

P P ~ P g„+l („-!)!• 


ln0 = ^2b n — 

71=1 ' L - 

b n = sum of all topologically 
distinct, labeled, doubly 
connected (irreducible) 
(n + l)-point graphs, subject 
to the labeling constraint 
that a vertex % must be 
connected to a vertex with 
label less than i. 



Table 1: Comparison of equilibrium Mayer expansion and random sequential ad- 
sorption. 



similar expansion for a lattice cooperative sequential adsorption model, while Tar- 
jus et al [0 obtained diagrammtic expansion from a very different method. 

The computation of the cluster integrals can be further simplified by noting that 
if a connected graph can be decomposed into several star graphs, the value of the 
connected graph is equal to the products of the values of star graphs. Because of 
this, we need only to evaluate the values of star integrals. An articulation point 
is a vertex belonging to the graph, such that when this vertex and its connecting 
edges are removed, the graph becomes disconnected. A star graph does not have 
an articulation point. The decomposed graphs of a connected graph are two or 
more graphs with the articulation point duplicated in each of the separated graphs 
when the articulation point is removed, with the edges still attached to the graphs. 
Figure ^ demonstrates some of the concepts. 

In the work of Dickman et al the series were obtained this way, using the 
existing cluster integrations in the fluid Mayer expansion. For the deposition of discs, 
the series was obtained up to S(4). For the squares of fixed orientation, the series 
is based on the result of one-dimensional cluster integrals 0] which were already 
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Figure 4: Some concepts of graph theory: (a) The labeling of the vertices consistent 
with RSA constraints; (b) this graph appears in equilibrium theory but not in RSA; 
(c) a connected graph with an articulation point, labeled A; (d) the connected graph 
(c) is decomposed into two star graphs (irreducible graphs). 



computed as early as in 1962 up to order 6. Note that for oriented hypercubes in d 
dimensions, a cluster integral is equal to the one-dimensional cluster integral raised 
to the power d. This fact was used to obtain two-dimensional results. Bonnier et 
al obtained results for d = 2, 3, and 4 for RSA of oriented squares, cubes, 
and hypercubes, using an extrapolation procedure from a lattice to continuum. In 
Table ^ we collect some of the series coefficients. 

The computational complexity of the RSA on continuum is perhaps more than 
(n!) 3 . One algorithm that we have implemented in this work consists of two steps. 
The first step reduces the labeled graphs to unlabeled graphs, directly from defining 
Eq. (|22"1). This gives us the RSA topologically distinct labeling factor. This step 
is purely graph-theoretic, and is independent of the detail of the problem. There 
are nfc=i(2 fc — 1) = 0(2" I 2 ) labeled graphs at order n. To reduce the graphs to 
unlabeled graphs, we consider all the n factorial possible permutations of the labels, 
taking the "largest one" as the canonical representation of the (unlabeled) graph. 
The largest one is the one with incidence matrix viewed as a big binary integer giving 
biggest value. Since there are 0(2™ I 2 ) graphs, each needs n\ transformations, the 
overall time complexity is 0(2 n / 2 n!). The above step reduces the number of graphs 
by a factor of about 1/n!. Much fast algorithms exist which generate each unlabeled 
graph exactly once. We used McKay's geng program j|7|, which greatly speeds up 
this part of the computation. 

The second step is to compute the integrals, computation of the cluster inte- 
grals can be done exactly for the discs only for n < 3. For the one-dimensional 
segments, it can be done in time of 0(n 4 n!). The result raised to the d-th power 
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771322713104711008093 
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24860884250598911650645 
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835568036857675195155997 
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29227671255970830546587445 



Table 2: The series expansion coefficient S{n) for the hard discs oriented squares 



[this work], and nearest neighbor exclusion on square lattice [36 . 



gives the (i-dimensional hypercubic result. For a given cluster integral, the set of fij 
gives a set of constraints. The integral is in fact equal to the volume of a convex 
polygon in an n-dimensional space. This n-dimensional space is partitioned into n\ 
regions, characterized by < x i2 < ■ ■ ■ x in where (i l3 i 2 , ■ • ■ , i n ) is a permutation of 
(1, 2, • • • , n). Once this extra constraint is given, the integration limits can be writ- 
ten down explicitly The integrals in each sector can be evaluated symbolically 
in polynomial time. 

This lastest program gives only two extra new terms after 35 hours of CPU time 
on a 600MHz Alphastation. The Dickman et al [Q and Bonnier et al [Q results 
were reproduced in 5 seconds. To obtain the next term requires about 10 3 times 
more in computer time, which is not possible within availability of our computer 
resources. 
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2.5 Other formulations of series expansions 



Beside the rate equation approach described in the previous three subsections, there 



are a number of other methods to obtain the series. Baram and Kutasov |32| and 
Baram and Fixman |55[] used Ising spin notation and master equation to relate 
expansion coefficients to a counting problem of lattice animals as follows. The 
particle density is expanded as 

PW = E(-!) n+ V^ -r^-. (26) 

n=l 

For the nearest-neighbor exclusion model, the coefficient a n is the number of con- 
nected lattice animals with n points that can grow from a single point. Unlike the 
rules embedded in the Eq. (|16]) for S(n) the new lattice animals are generated from 
perimeter sites only. 

We can get this result from the rate equation point of view as follows. Let 
P(W) be the probability of a set of connected sites W for each of which a particle 
can be put in for sure. A large set G must be known to be empty. However, the 
correspondence between W and G is many to one. The rate equation for P(W) is 
then 

dP{W) 



dt 



-m 



P(W)-J2P(WU{x}), (27) 



where m = \W\ is the number of sites in W, and x runs over the perimeter sites of 
W (for the nearest neighbor exclusion model). Let us transform t to y — 1 — e - *, 
and Q(W, y) — (1 — y)~ m P(W, t), then the new equation is 

iQ(W) -- Y. WV'UM). (28) 



dy 



perimeter of W 



This naturally gives Eq. (^) and rule for a n . 

Dickman et al [33| gave an operator formulation of the master equation. The 



formal solution of the master equation is used to derive expansion rules identical 
to the rate equation approach. For the continuum problem, Schaaf and Talbot | 29| 



gave small density expansions based on results in scaled particle theory P5 |. The 
final result is equivalent to the time expansion, but the procedure is quite different. 

Consider the average fraction <fi of surface available to the center of a new disc, 
(f) = dp/dt. One has 

<\> = 1 - S 1 + S 2 - S 3 + 5 4 , (29) 

with 

Sn = nJ.I I ■ J dxxdx 2 ■ ■ -dx n A n (x 1 , ■ ■ ■ ,x n ) p (n) (^i," • • ,x n ), (30) 

where A n (xi, ■ ■ • , x n ) is the area common to the exclusion circle of n particles ad- 
sorbed at the positions defined by the vectors x\, x 2 , ■ ■ ■ , x n \ i.e., it is the area of 
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Figure 5: Series results for the RSA coverage on a square lattice of the monomer 
with nearest neighbor exclusion model. The curves labeled 19, 20, and 21 are the 
direct time series truncated to 19, 20, and 21 orders, respectively. The line labeled 
"[10,10] Pade" is N = D = 10 Pade approximant in variable y, defined by Eq. (j32j), 
with b = 1.05. For convenience of viewing the whole range of t, we plot the curves 
with variable 1 — exp(— t). 



the set D(x\) PI D(x 2 ) • • - fl D(x n ). The multi-variable function p^ n \x 1 , ■ ■ ■ ,x 2 ) is 
the general n-particle distribution function, which is proportional to the probability 
that a particle will be found in the elements of surface dx\ at xi, a second in dx<i at 
X2, etc. Schaaf et al observed that the leading order in a density expansion for S n is 
equal to a corresponding equilibrium system. This allows them to find exact expres- 
sion for S n correct to third order in density [^9[ |49], [50], |51|| . Tarjus et al [§|| derived 



Kirkwood-Salsburg-like equations for p^(xi, ■ ■ ■ ,x n ). From these equations, a di- 
agrammatic expansion in density is obtained. Some what opposite in spirit to the 
rate equation approach, Given |Q formulates nonequilibrium problems as quenched, 
multi-species equilibrium problems. He also gave the topological reduction in graph 
expansion, showed clearly the parallel between RSA and equilibrium system. Caser 
and Hilhorst presented series expansion directly for the jamming coverage. The 
approach is interesting, but the method is not very accurate. 
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2.6 Analysis of the series 



The series obtained for the coverage in t has a finite radius of convergence, see 
Fig. [|. The fundamental question is whether we can give good approximation valid 
for the whole time domain t G [0, oo). Although we do not have rigorous proof, the 
numerical procedure which we shall discuss below indicates a positive answer. 

The standard technique of extending the radius of convergence is the Pade ap- 
proximation method |22], |24]]. Given a series f(x) to order L, we determine two 
polynomials Pn(x) and Qd(%) of degree N and D respectively, such that 

m-^rr\ = 0(x L+1 ), N + D < L. (31) 
Qd{x) 

It is a powerful way of extending the domain of convergence of the original series. 
The polynomials P/v(x) and Qd(x) with Qd(x) = l + b\X+- ■ ■ is usually determined 
uniquely for the given f(x). 

To accelerate the convergence further, new variables are introduced, e.g., 

y=l-exp(-6(l-e - *)). (32) 

The functional form is chosen in such a way so that the series in the new variable y 
most closely resembles the asymptotic behavior of the coverage 9(f) at large t. For 
small t, y is proportional to t. The above specific form is encouraged by the exact 
solution of the one-dimensional dimer problem. In fact, we have 9 = y with 6 = 2 
in such case (see Eq. (|J)). The convergence among various Pade approximants is 
improved greatly by the transformation. This form of transformation works well for 
all lattice models that we have studied. For continuum models, Dickman et al p3| 
considered 

v = 1 ~W+u (33) 

for the RSA of discs in two dimensions and 

y 1 + bt K J 

for the oriented squares, which match the asymptotic law of t~ l l 2 (known as Feder's 
law 0) and log t/t law (Swendsen ||), for discs and oriented squares, respectively. 

Other methods of analysis were also used. For example, Fan and Percus |53| 



used exact results of solvable systems as reference systems which in some sense are 
close to the actual system to speed up the convergence. The coverage of the nearest 
neighbor exclusion model on a square cactus fractal lattice is |)3 



0W)) = \ + \*-\&-*)\ (35) 
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where variable z is related to the original variable t by 

dz' 

/o i + icT-z 



1 - e~* = 3 /_ - , 0/1 , 3 . (36) 



We note that the coverage 9 is a truncated power series in z for the cactus system. 
The idea is that if we use the same variable z, relating t through Eq. (0), the 
nearest neighbor exclusion model on square lattice should have a series in z with 
much improved convergence. 

Baram et al j32, ^| used Euler transform 



1 — u , , 

u = e~\ cwl, (37) 



c + 1 — u ' 



and Levin convergence acceleration method [|54| for the series analysis. It appears 



that the quality of the results is comparable to other methods. 

The algebraic manipulations can be done with great ease by symbolic packages 
such as Maple or Mathematica. Mathematica contains function to find the Pade 
approximation of a series. The results of the analysis can then be plotted within the 
software. Here is an example of simple Mathematica code to form the series, do the 
transformation, call the Pade function, and finally evaluate the jamming coverage. 
The SetPrecision[l .35,40] command makes the subsequent computation in 40- 
digit precision. 

«Calculus'Pade' 
maxord=17; 
P[0] = 1; 
P[l] = -4; 
P[2] = 28; 
P[3] = -268; 

. . . (omitted) 
P[17] = -6058617368871081964076; 

theta = 1 - Sum[P [i] *t~i/i ! , {i ,0 ,maxord}] + [t] ~ (maxord+1) ; 
b = SetPrecision[l. 35,40] ; 

f = theta /. t -> - Log[l + Log[l-y]/b] + [y] " (maxord+2) ; 

yinf = 1 - Exp [-b] ; 

pd = Pade[f , {y, 0, 8, 8}] ; 

thetainf = pd / . y -> yinf 

The Pade approximant for the coverage of the dimer on square lattice is given 
in the pd variable as 

(2.962963 y + 0.03206897?/ 2 - 2.195246 y 3 - 1.073721 y 4 + 0.9207869 y 5 + 0.5556586 
-0.04386743 y 7 - 0.05303456 y 8 ) / (l + 1.733045 y - 0.2568919 y 2 - 1.942572 y 3 
-0.5852424 y 4 + 0.7908992 y 5 + 0.4421557 y 6 - 0.0493306 y 7 - 0.0513337 y 8 ) , (3 

where y is given by Eq. (|32|) with b = 1.35. The result is stable against variation 
in b. Error can be estimated from the convergence of various Pade approximants. 
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0.3641327 
0.3641326 
0.3641325 
0.3641324 



0.3641323 
0.3641322 
0.3641321 
0.3641320 

0.0 0.5 1.0 1.5 2.0 

b 

Figure 6: Pade approximant estimates for the jamming coverage 9(oo) as a function 
of the transformation parameter b, for the monomer RSA with nearest neighbor 
exclusion on a square lattice. The transformation is given by Eq. (0). (from 
ref. [§§). 




The above Pade approximant is accurate to 10 5 for all t > 0. Such a high degree 
of accuracy embedded in a simple formula is remarkable. 

We set b to be a free variable so that when different orders of Pade approximants 
are applied to 0(y), a crossing region between different orders of Pade approximants 
is to be located so as to give a good estimate of jamming coverage 9(oo). Some 
of the new series results and analyses are presented in Gan and Wang [[36|]. As an 
example, Fig. |6| shows the crossing region of Pade approximants of orders [N, D] 
with 19<iV + .D<21,.D>8,iV>8, for the series of monomer RSA with nearest 
neighbor exclusion on a square lattice, giving an estimate of 9(oo) = 0.3641323(1), 
where the last digit in parentheses denotes the uncertainty in the preceding digit. 
This estimate is in good agreement with the estimate of 6(oo) = 0.3641330(5) by 
Baram and Fixman ||35||. Analyzing the series using the square cactus as a reference 



model |j5"3"H , through the transformation Eq. (|36|), we obtain 9(oo) = 0.364132(1), 
less accurate than that using transformation Eq. (^). All these estimates from the 
series analysis agree well with the simulation result of 6(oo) = 0.36413(1) [p5| . 

The estimates for the jamming coverage on continuum systems are much less 
accurate, due to limited number of terms available in the series. In Fig. [7] we present 
the Pade estimates for the oriented squares as functions of the parameter b. Best 
convergence occurs around b = 1.3, given 9(oo) = 0.5623(4). The additional new 
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0.550 I— 
1.20 



1.25 



1.30 

b 



1.35 



Figure 7: Pade approximant estimates for the jamming coverage 9(oo) as a function 
of the transformation parameter b, for the oriented square. The number [N, D] is 
the order of the Pade approximant. The transformation is given by Eq. (|34|). 

terms do not improve the result of Dickman et al very much. 

In Table we collect the results of jamming coverage compared with Monte 
Carlo results. It is clear, if the jamming coverage can be estimated with great 
accuracy, so is the whole function 9{t). The results for continuum models are less 
accurate. Nevertheless, we still obtain quite reasonable estimates. 



3 Monte Carlo simulation 



A large variety of RSA models have been simulated. Lattice models are easy to 
simulate with good accuracy. Random dimer filling is given by Nord and Evans p6| 
and Wang and Pandey [Q. Nearest neighbor exclusion models are simulated by 
Meakin et al Square blocks of n x n are simulated by Nakamura ||57|| , and by 
Barker and Grimson [p£| . Square blocks with an emphasis of approaching to the 
continuum limit are studied by Privman et al [p9|] . 

On continuum, hard-disc system is simulated by Akeda and Hori 



50 1, Finegold 



and Donnell HI, Tanemura 61 1, Feder 101, Hinrichsen et al 62 , Meakin and Jullien 



63], and Wang [S4|. The oriented squares were of interest [BC] because of Palasti 



conjecture 



which says that the jamming coverage for the two-dimensional ori- 
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Mnrlpl 
lVl ULLC-L 


Series 


Monte Carlo 
0(oo) 




Pade 6(oo) 


NN (square lattice) 
dimer (square lattice) 
NN (honeycomb) 
dimer (honeycomb) 
NNN (square lattice) 


21 
18 
24 
22 
14 


0.3641323(1) [36 
0.906823(2) 36 
0.37913944(1) [36 
0.8789329(1) p> 
0.186985(2) [35 




0.36413(1) ! 

0.906820(2) 

0.38(1) 

0.87889 

0.186983(3) 


55 
56 
28 
PB] 
59 




hard discs 
oriented squares 


5 
9 


0.5479 |3| 
0.5623(4) [this work] 


0.5470690(7) [H 
0.562009(4) [£6J 



Table 3: Some of the most accurate jamming coverages by series analysis and Monte 
Carlo computer simulation. NN stands for nearest neighbor exclusion, NNN for 
nearest and next neighbor exclusion. The order n is the highest known term in the 
time expansion of 0(f). The numbers in square brackets are references. 



ented square is equal to the square of the one-dimensional line- segment jamming 
This conjecture turns out to be false 



66 



coverage. 

Other more elaborated geometries are considered such as randomly oriented el- 
lipses and lines |68|, |6£ 



randomly oriented rectangles [70 



and other anisotropic 
particles [J7T|] , spheroids |72[ , mixtures of different sizes |53], [73|, [73J] , The interests are 
the jamming coverage and the asymptotic laws for large time. Clearly, these com- 



plicated problems are difficult for systematic series expansions |5l| . Lattice models 
with dimers on square lattice deposited with different probabilities in two orienta- 
with line segments f76j, mixtures |77], |78|, ff9j , and random walks |56|, |8C] 



tions 75 



are also studied. 

We mention here briefly two related generalizations of standard RSA. The first 
of these is multilayer RSA ]8T| . The jamming densities as a function of the height 
show interesting power-law behavior 



s2 



diffusional relaxation HO], B3L |8~4l, ISo" 



The second generalization is RSA with 
, BSl. When diffusional relaxations are in- 



troduced, many of the intrinsic properties of RSA changes. Most of these generalized 
models are studied by Monte Carlo computer simulation. 



3.1 Simple algorithms 

Monte Carlo simulation method has the advantage of very little programming in- 
vestment, quick and robust results. Random sequential adsorption, especially lattice 
models, can be simulated on a computer rather easily. The basic variables (declared 
as an array) are the occupations of the lattice sites. At each time step, a site x is 
picked at random, if some neighborhood D(x) is empty, we put down the particle by 
reset the array values in these entries of involved sites. Note that one is simulating a 
time-discretized version of the RSA model, while the rate equation approach consid- 
ers continuous time. The discretization step is 1/L d where L is the linear dimension 
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of the system. However, discretization effect is unimportant since we can simulate 
very large systems easily. 

With continuum models, since we can not have occupation variables like the 
lattice models, more tricks are needed for efficiency. The coordinates of the particles 
are stored to register the locations. To speed up the checking for overlapping, 
"coarse-grained occupation variables" are needed to indicate where are the particles 
in a given range, a technique borrowed from molecular dynamics simulation. For 
example, in RSA of discs, the total L x L surface area is partitioned into regions of 
squares of size 1 x 1. A disc at (x,y) belongs to the unit square at ([x\, [y\), where 
the notation \_x\ represents the floor or integer part of x. An array of linked lists of 
disc coordinates is used to represent the deposited discs. With this data structure, 
it is suffice just to check the discs located at the center and eight surround squares 
for overlaps. 

The above simple methods have a severe shortcoming that the study of late-stage 
process is very time-consuming. When the jamming configuration is about to be 
reached, much of the area is blocked. Only tiny disconnected regions are available for 
further deposition. Since we choose the sites at random with a uniform probability 
distribution, the available sites are hard to find. More advanced methods cleverly 
discover these regions and suggest deposition only in regions which are potentially 
possible for adsorption. 

3.2 Event driven algorithms 

First we consider the lattice models. On a lattice, since there are only finite number 
(of order N = L d ) of possibilities for deposition, each with equal probability in 
the original definition of the model, we can classify the deposition possibilities as 
possible and impossible attempts. For example, for the nearest neighbor exclusion 
model, each site is associated with a deposition attempt. The possible sites are 
those with center and four neighbors empty; the impossible sites are those with any 
one of the five sites occupied. Let us assume that we have N a possibilities for the 
available type and iVj for the impossible type. Clearly, N a + iV; = N. In the normal 
simulation, we hit the first available type with probability p = N a /N and other type 
with 1 — p = Ni/N. Consider a total of k attempts such that the first k — 1 failed 
to deposit a particle because one hits already blocked sites and the A;-th attempt 
is successful and is deposited. What is the probability that such an event occurs? 
Each failed attempt occurs with probability (1 — p), so k — 1 consecutive fails have 
probability (1 —p) k ~ 1 - Since each failed attempt does not change the current state 
of the system, the next attempt is totally independent of the previous attempts. 
The successful attempt has probability p so the total probability is 

P k =p(l-p) k -\ fc = l,2,.... (39) 
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Imagine that we do not actually do the k — 1 failed attempts, and pick one at random 
from the available sites and deposit the particle. From the above discussion, this 
is equivalent to k attempts in units of the original time scale. However, k is not a 
definite number, rather it follows the probability distribution, Eq. (|39|). 

Event-driven type of algorithms was invented long-time ago for the simulation 
of Ising models |89| , which is known as N-fold way. An event-driven RSA algorithm 
64j , |66l on lattice is as follows: pick deposition site from the list of available sites 



at random. Increase the time t by k/L d (L d attempts are usually defined as unit 
time). Update the configuration and the list. 

The random integer k can be generated by transformation method as follows: 



1 + 



ln£ 



ln(l — p) 



(40) 



where £ is a uniformly distributed random number between and 1 . Some program- 
ming efforts are needed to update the list efficiently. After one deposition, not only 
the current site but also some other sites in the neighborhood of the current site are 
no longer available. This should be correctly reflected in the list. The deposition 
also changes the value p. 

Perhaps the best application of the technique to lattice model is the RSA of 
polymer chains [Q. Due to a large number of conformational states of polymer 
chains, the deposition becomes very slow in the late stage. An event-driven method 
is, therefore, essential to obtain long-time results. This is accomplished by identi- 
fying the early part of growing chains, i.e., the partial chains and then classifying 
them periodically as "available" or "not available" for the deposition in an itera- 
tive fashion. The subsequent depositions are then made starting from the available 
partial chains choosing at random. The list of the partial chains is limited by the 
available computer memory and governs the speed of the program. It appears that 
there is an optimal list size for a given length of a chain. 

Highly accurate Monte Carlo simulation results were obtained for oriented squares 
|66| and for discs |M] on two-dimensional continuum. In the RSA simulation of discs, 



we incorporated two important ideas: (1) divide the surface into small squares and 
make deposition attempts only on squares that are not completely blocked; (2) sys- 
tematically reduce the sizes of the small squares after some number of attempts and 
re-evaluate the availability of the squares. 

In this more elaborate method, we make attempts only on squares which are 
potentially possible for depositions. Thus the core part is an algorithm which iden- 
tifies correctly whether a square is available for deposition. Since the diameter of the 
discs is a = 1, each disc excludes an area of circle of radius a for further deposition. 
If some area is completely covered by a union of the exclusion zones of discs, that 
area is not available for deposition. We begin by classifying all the squares of a x a 
with a = 1 as available or unavailable squares. The available squares are put on a 



23 




Figure 8: The six different situations for discs covering or not covering a square. 

list. Deposition trials are taken only on the squares in the list, i.e., a square in the 
list is chosen at random, deposition is attempted with a uniform probability over 
the square. After certain number of trials (we used 5 x 10 5 ), we rebuild the list, now 
with squares of size (a/2) x (a/2), checking only those squares on the old list. Each 
old square is subdivided into 4 smaller squares. This shrinking of the basic squares 
helps to locate even extremely small available regions. They will be hit with very 
small probability if the area is always lxl. Because of integer overflow, the size 
of the squares is not allowed to go arbitrarily small, but the shrinking is stopped at 
a = 2~ 15 . Thus the smallest square has an area of 2~ 30 « 10~ 9 . 

For the standard RSA algorithm, each trial increases the time t by I /A, where 
A = L 2 is the area of the total surface. In the even-driven algorithm, where only 
the potentially successful depositions are tried, the clock must tick faster in order to 
compensate for not making attempts in the completely unsuccessful area. The time 
increment for each trial on the available squares is given by k/A with k generated 
using Eq. (]4"0D, with p being the ratio of the area on which we try our deposition to 
the total area; it is equal to a 2 /L 2 times the number of available squares. 

An important piece of code of our RSA simulation program is the identification 
of the squares which are fully covered by the excluded area of discs. For these 
squares we are sure that depositions on them will be unsuccessful, and thus they 
will not be on the list of candidates for trials. 

The algorithm that we have devised is as follows. Written as a C programming 
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language function, it returns a nonzero value if the square is fully covered and a 
value otherwise. Part (1) to (5) is executed in the order given. 



1. Find relevant discs to the current square. A disc is relevant if the square 
overlaps with the excluded region of the disc. Note that a disc has a diameter 
a and its excluded region is a concentric circle of radius a. 

2. If all of the four vertices of the square are inside a single-disc excluded region, 
then it is already fully covered (Fig. ||a). Function returns with value 1. 

3. For a full coverage, each vertex of the square must be at least in one of the 
excluded region of the relevant discs. If any one of the vertices is not covered 
at all by any disc, then the square is not fully covered (Fig. ||b). Function 
returns with value 0. 

4. At this point, at least two discs are involved if the function is not returned. 
Mapping the edges of the square to a one-dimensional line between and 4, 
we find out all the line segments which are covered by the excluded region of 
discs. If the four edges of the square are completely covered by the discs, and 
we have exactly two discs, then the square is fully covered (Fig. ||c), function 
returns with value 2. If there are segments not covered by discs, then we know 
that the square is not fully covered (Fig. |8|d). Function returns with value 0. 

5. If the function is not returned, we know that there are at least three discs. 
Even if all the edges are covered, the square can still contain uncovered area if 
three or more discs are involved. Three or more discs can form holes. To make 
the last check, the coordinates of the intersection points of circles (with radius 
a) are calculated. In order for the interior of the square to be covered, all the 
intersection points of any pair of discs which lie inside or on the square edges 
must be covered by a third disc. Return with the number of relevant discs if 
the above condition is satisfied (Fig. |S]e); return value if not (Fig. ^|f). 



The order of various checks is arranged in such a way so that the most frequent 
and easy situations are checked first. Even though the worse case computational 
complexity goes as n 2 , where n is the number of relevant discs, a definite conclusion 



can be made typically much earlier. Using this algorithm, Wang |64| obtained 
accurate jamming coverage 6{oo) = 0.5470690(7) and confirmed to a high accuracy 
the Feder's asymptotic law 6{t) fa 9(oo) — ct^ 1 ^ 2 . 



4 Conclusion 



We introduced the method of series expansions for both lattice models and contin- 
uum models. The computation of the series for the lattice model can be reduced to 
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a counting problem. An efficient implementation requires full use of the symmetry 
of the problem. For continuum models, systematic expansion in terms of a class 
of graphs are available. Both the problem of reducing labeled graphs to unlabeled 
graphs and of computing the cluster integrals are computationally difficult problems. 
Thus, the series for continuum systems are rather short. The series can be analyzed 
through variable transformation and Pade approximation. Monte Carlo simulation 
method is versatile, and implementation is straightforward. However, efficient algo- 
rithms are available with more sophisticated data structures. The challenge to the 
series expansion and efficient Monte Carlo simulation for RSA is to apply to more 
complicated models which may describe better real experiment situations. 
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